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ABSTRACT 


In the numerical integration of orbits by a multistep 
process, it has been suggested that the application of a local 
error control may increase the efficiency of the integration 
without any significant loss of accuracy. This study de- 
velops methods of controlling the local error and examines 
the effects of these controls on the integration. It is shown 
that for various orbit types, controlling the local error can 
optimize the integration while maintaining the desired 
accuracy. 
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LOCAL ERROR CONTROL AND ITS EFFECTS 
ON THE OPTIMIZATION OF ORBITAL INTEGRATION 

by 

C. E. Velez 

Goddard Space Flight Center 


INTRODUCTION 

Many computing systems used to determine the orbits of artificial earth satellites require nu- 
merical integration. Due to recent advances in the theory of perturbations and the concomitant in- 
crease in the complexity of the mathematical model, highly efficient integration techniques are 
required. Such a technique is described in this paper, based on the concept of ,f local error control.” 

Consider the orbital equations 


x = W 1- + P “ x ° ’ (*) 

in three space variables, where lx| = (x 2 +y 2 + z 2 ) 172 , P is the perturbation function, and ^ is a con- 
stant. We assume throughout thatP is fairly complicated, so that the efficiency of the integration 
is proportional to the total number of derivative evaluations. Also, for simplicity, we assume that 
p is independent of the first derivative x, i.e., p = P(t, x). 

The numerical method under consideration is a multistep process, i.e., a method that approxi- 
mates the solution of Equation 1 by using formulas of the form 



n - k , k + 1 , * * • , 


( 2 ) 


where h is the stepsize, x. = x(t 0 + ih) , and a., /?. are constants. Formulas of this type define a 
linear k-step method (Reference 1); they are usually derived by integrating a polynomial approxi- 
mation of the derivative x . 

Associated with Equation 2 is a local truncation error of the form 

R n = R(t n ) = Ch p+2 x<p +2 > (O , (3 ) 
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where <f is contained in [t n _ k# t n ] , c is a constant, and pis an integer called the "order" of the 
method, all depending on k, the number of "backpoints" used. 

Soar has suggested (Reference 2) controlling the local error by varying the order and stepsize 
while integrating Equation 1 with this object: It is clear (for example) that the magnitude of R n is 
sensitive to variations in p or h. Suppose that during the integration the magnitude of R n becomes 
insignificant relative to the required accuracy of the calculations being performed. Then increasing 
h or decreasing p, separately or in combination, may increase integration efficiency without sacri- 
ficing accuracy. The purpose of this study, then, is to develop automatic methods of controlling 
the magnitude of R n during integration by varying the parameters p andh, and to examine the effects* 
of these controls on efficiency and accuracy. Before discussing methods of estimating and con- 
trolling the local error, we formulate a commonly used integration model that was used to obtain 
the conclusions given under "Numerical Results." 


THE INTEGRATION FORMULAS 

In Equation 2, if /3 0 ^ 0, then knowledge of the solution x n is required on both sides of the 
equation and cannot, in general, be explicitly solved. However, equations of this type (closed form) 
have smaller associated truncation errors as well as desirable stabilizing characteristics (Ref- 
erence 3). The well-known predictor-corrector algorithm uses formulas of this type by first com- 
puting an initial (predicted) approximation of the solution using a similar formula with /3 0 = 0, then 
using the closed form iteratively until convergence is achieved. It can be shown that for sufficiently 
small h, the successive corrected values obtained by this process converge to the unique solution 
of the closed-form equation, provided that the function being integrated is sufficiently smooth. 

Consider now the predictor-corrector formulas that are derivable from Newton's backward 
difference interpolation polynomial (Reference 1, pp. 291-293— see also Appendix A): 


Predictor: x n + 2x . + x 9 - V 2 x 

n n 1 n J. n 


Corrector: x^ - 2x . + x - V 2 x 

n n 1 n ^ n 


h 2 [i +^- V 2 V 3 V 4 V s + ••• +(T q v q ] x n -i , (4) 

h2 [ 1 _ V + ^ V 2. 5 1 5 V 4__1_ V 5 + ... + a* V*)] £ n . (5) 


These are the Stormer-Cowell formulas expressed in terms of backward differences. The 
local truncation errors associated with the formulas are given by Equation 3, where P = q + 1 and 
C = crq+1 or (j q * +1 . If the parameter q is fixed (and hence the order p), then Equations 4 and 5 can be 
written in the form of Equation 2 if the backward differences are expressed in terms of the 
ordinates x r 


* Local error control does not generally yield a quantitative measure of the accumulated error. However, a qualitative control is pos- 
sible and in most cases sufficient. For a discussion on accumulated error estimations, see Reference 1. 
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In particular, for q = 5 there are the Stormer-Cowell sixth-order formulas: 


h 2 


v2 *n = 240 [^17 x n _j - 266 x„_ 2 + 374 x n _ 3 - 276 x n _ 4 + 109 x n _ s - 3 x n _ g ] 


V2x n = 240 [ 18x n +209 x n-l +4x n-2 + 14X n-3- 6x n-4 + X n-s] 


(4) ' 

( 5 ) ’ 


The coefficients in Equations 4 1 and 5 T clearly depend on the choice of order, so that varying 
the order during integration would mean producing a new set of coefficients for each p. Therefore, 
formulas such as Equations 4 and 5 (where variations of the order can be made simply by varying 
the number of terms retained) were used in this study.* Before any of the above formulas could be 
used, a set of ’’starting” values of the solution must be computed^ For example, for Equations 4 T 
and 5 T , if the values x. (and hence x^, i = 0, 1, • * • , 5 are known, then Equation 4 f could be used 
to obtain a ’’predicted” value x p (n = 6) and Equation 5 f to obtain the successive corrections x 6 C ’, 
j = 1, 2, 3, • • • , until convergence to the same criterion 5 is achieved, i.e., until 


The process could then be repeated with n = 7, 8, • • ■ . 


LOCAL ERROR ESTIMATION 

Any considerations concerning the control of the local error depend on the capability of ob- 
taining a reasonable estimate of Equation 3 during the integration. A widely used approximation is 
based on the fact that if f(x) is an n-times continuously differentiable function, then there exists a 
quantity <9, 0 < 8 < 1, such that 


A n f ( x ) 
(Ax) n 


d n f 
dx n 


[x + nd (Ax) J , 


i.e., if Ax is sufficiently small, we can approximate a high-order derivative by a high-order dif- 
ference. In particular, we can write 


R n = R ( t n) = ChP +2 X<P +2 >(£) = Ch P + 2 X< P >(0 * Ch 2 V p X n - (6) 


*In actual practice, a modification of Equations 4 and 5 known as the "summed” form of the integration formulas was used. This modi- 
fication is formulated in Appendix A. 

^Experimentation indicates that recently developed high-order Runge-Kutta type formulas are particularly suited for such purposes (Ref- 
erence 4). 
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If suitable predictor-corrector formulas are used, another technique is available: comparing 
the predicted value x n p with the finally accepted "corrected" value x n c . In particular, it can be 
proven (Reference 1) that, for Equations 4 and 5, 

E n = K(x n C - X P) , ( 7 ) 

where K is a constant. This technique is known as Milne T s method of estimating the local error. 

LOCAL ERROR CONTROL 

Consider now a specific technique which makes it possible to control the magnitude of the local 
error during integration by varying parameters p and h. Let T x and T 2 be tolerances specifying 
upper and lower bounds on the local error so that, for any n, the local error must satisfy 

T 2 - K| - T 1 ' where T 2 < T, . ( 8 ) 

Controlling the local error by varying the order can then be accomplished by determining 
whether this condition is satisfied for each value of n; if for some n, R n <T 2 then decrease the order; 
if R n >T 1? then increase the order. However, varying the order alone is generally not enough for ef- 
fective control. For example, if, for some n, R n >T and h is large, it may not be possible to satisfy 
Equation 8 with any p; in which case the stepsize must be decreased. (Also, the danger of numer- 
ical instability increases considerably with large p; see References 5 and 6.) On the other hand, 

R n <T 2 and p small would indicate that a larger stepsize could be used. 

Controlling the local error by varying both stepsize and order can be accomplished as follows: 
Let L l and L 2 be limits specifying the desired upper and lower bounds on the order, i.e., 

L 2 < p < L x where L 2 < L x . (g) 

Then, if at some point during the integration, p >L X , decrease the stepsize; if p <L 2 , increase the 
stepsize. 

In a local-error control so designed, the parameters T x , T 2 , L x , L 2 completely govern the 
degree and type of control; i.e., T x and T 2 can be selected so as to effect any degree from no con- 
trol to continuous step-by-step control; likewise, l x and L 2 can be selected so that the control in- 
volves varying the order alone, varying the step alone, or varying both. (For example, if L x = l 2 , 
the control depends solely on stepsize variations.) 

Changes in stepsize magnitude during integration (unlike changes in the order) are usually hard 
to accomplish, since they need a "memory" of equally spaced points at every step during integra- 
tion. For this reason a common technique is "halving- doubling," where an increase or decrease in 
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stepsize is half or double the current stepsize. Then increasing the stepsize presents no problem 
if a sufficient number of backpoints are retained. Decreasing is usually accomplished by interpo- 
lation or by a single-step method. 


OPTIMIZATION OF STEPSIZE 

Restricting the variations of stepsize h by a constant factor does not generally yield "optimal” 
stepsizes;* hence the initial choice of interval could have a substantial effect on the total number 
of integration steps. For example, suppose the halving-doubling method is being used and the order 
is held fixed. If for some n, R n >T 1 , the stepsize must be decreased. Let h opt be the optimum 
(largest) value of h for which Equation 8 is satisfied. It might then happen that h/2 < h opt < h. But 
in this method, h/2 would be used, although this would necessitate more integration steps than if 
the optimum stepsize were used. Hence a variation of h that would better approximate its optimum 
value is desirable. One technique available (see Reference 7 for details and applications in the 
case of single-step methods) is to compute the stepsize using the local-error estimate, i.e., letting 
a be the "allowable” local error for each step, where T 2 <cr<T l . Variations in h can then be com- 
puted from the relationship between a and R n *= ch p+2 x< p+2 > (/) : 

r q ~|l/(P + 2) ^ 

h ° pt ' Lcx^ 2 ) (£)„ 

so that, if cr <R n , the stepsize is decreased; if R n <a the 
approximately optimal with respect to the choice of o-. 


crh p+2 

R 


l/(P+2) 


( 10 ) 


stepsize is increased; the variation being 


OPTIMIZATION OF ORDER 

Varying the order by a constant factor, like varying the stepsize, need not yield the "optimal" 
order. 1 " For example, suppose that a variable-order, constant-step control is being used and the 
order is varied by ±1. If, for some n, | R n | <T 2 ,the order must be decreased. Let p opt be the opti- 
mum (smallest) value of p for which Equation 8 is satisfied. Then the following situation could exist; 

Popt < P " 1 < P • 

In this method, p - 1 would be used, although this would result in more calculations than are 
required at the optimum order. Also, using a variable-step control in addition to varying the order 
may cause: 

Popt < L 2 - P - 1 , 


*The largest stepsize that allows a prescribed local error at a given point. 
^The smallest order that allows a prescribed local error. 
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in which case the stepsize would be increased if p opt were used. One method to obtain the optimum 
order is to test the local error for various orders until the optimum is established. For example, 
if the local error estimate used is given by Equation 6, then ch 2 V? x n can be tested for various p 
until the smallest value satisfying Equation 8 is found. 


GENERAL EFFICIENCY CONSIDERATIONS 

Yhis section discusses some possible effects of the above mentioned controls on the efficiency 
of any integration, which will generally depend on attaining the following objectives: 

A. Minimizing the number of integration steps, 

B. Minimizing the number of corrector iterations required at each step, 

C. Minimizing other computational efforts required by the integration formulas being used, 
such as retaining only the significant terms in Equations 4 or 5 during integration, 

D. Minimizing the effort and time involved in computations to control the local error, such as 
producing the required ’’memory” when changing the stepsize. 

Three types of control are now considered. 

1. A variable-order, fixed-step control: 

Since varying the order involves only varying the number of terms retained, this control has 
little effect on attaining A. Equation 7 however, expressing the relationship between the local 
error and the ’’predicted-corrected” difference, indicates that controlling the local error may sub- 
stantially affect the total number of corrector iterations. In particular, an increase in order at 
some point during integration where the local error is increasing could minimize any increase in 
the number of required iterations. This control also has an obvious effect on attaining C, and also 
that given a stepsize, the order required to satisfy a given criterion on the local error is automat- 
ically selected. 

2. A fixed-order, variable-step control: 

This control affects the number of corrector iterations for the same reason as control 1. It 
evidently affects the attainment of A or D, so one of them must be sacrificed. The gain in efficiency 
from larger stepsizes must be weighed against the cost in changing the stepsize. Furthermore, 
there are the two methods for varying the stepsize: halving- doubling, and step computation using 
the local-error estimate. In both cases the necessary ’’memory” could be obtained by interpolation, 
but in the second method all the backpoints must be computed when decreasing and when increasing 
the stepsize, as compared with opposed half, or none, for the first method. However, the second 
method may be the more effective of the two in minimizing the number of integration steps. This 
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control has the further advantage that, given an order, it automatically selects the stepsize (possibly 
optimum) required to satisfy a given local error criterion. 

3. A variable-order, variable-step control: 

The control combines the effects of controls 1 and 2. It has one possible disadvantage: in 
allowing the order to vary before changing the step, it may cause smaller stepsizes, thus opposing 
objective A. But this control has one clear advantage: it automatically selects both the stepsize 
and order required to satisfy a given local error criterion. 


EFFECTIVE ERROR CONTROL IN ORBITAL INTEGRATION 

In general, the effectiveness of a local error control during the integration of a particular orbit 
depends on the degree and rate of change of the derivative x (p+2) (£) (and hence the local error), 
during a revolution for a fixed p and h. This change is governed by the orbital parameters a (the 
semi-major axis) and e (the eccentricity) (Reference 2). In particular, for orbits with eccentricity 
near zero, the local error variation will probably be small over a revolution, and as the eccentricity 
becomes bounded away from zero, this variation becomes more pronounced. Moreover, the larger 
the semi-major axis for a particular satellite, the slower the rate of change of the derivative. 

Suppose that during the computation of a particular orbit, some local error criterion is to be 
satisfied over the entire range of integration, which may involve many revolutions. That is, sup- 
pose that the local error is to be bounded from above, so as to restrict (at least qualitatively) the 
propagation of truncation error. Suppose also that the integration method used is of some fixed 
order (as is generally the case when formulas such as Equations 4 T and 5 f are used), and that the 
stepsize is to be specified. The error criterion could then be satisfied, without serious loss of ef- 
ficiency, by integration over a single revolution with various stepsizes— choosing the largest step- 
size that satisfies the criterion over the whole revolution. For orbits with e ^ 0, this process 
could yield the optimal stepsize for the given order. On the other hand, for orbits with e > e > 0, 
the stepsize selected in this manner would mean bounding the local error at its maximum value 
(e.g., at perigee), and could cause needless computation where the local error is smaller or at its 
minimum (e.g., near apogee). 

Consider the effects of a continuous variable-step control during integration. In both cases 
cited above (i.e., e = 0, or e > e > 0) an initial stepsize (possibly optimum) satisfying the local er- 
ror criterion would be computed automatically; furthermore, the stepsize would vary according to 
variations in the derivative. The numerical results show that in both cases, under suitable condi- 
tions, such a control has a considerable effect on the efficiency. In particular, if a sufficiently large 
range (T^ t 2 ) is used in Equation 8 (so that the local error is allowed to range between these two 
limits before any interval modification is performed), the gains in efficiency due to the larger step- 
sizes will overshadow any loss due to changing the stepsize. 

If, in addition to the above, the order to be used could be specified (as it generally can when 
formulas such as Equations 4 or 5 are used), then it may be possible to obtain an optimal stepsize 


I 



(or stepsizes)— order combination by integrating over a single revolution with various orders and 
letting the stepsize vary during the revolution. In particular, we could find that order which— 
together with a continuous step modification— gives the most efficient integration. Examples of 
such a procedure will be given in the numerical results. 

Finally, consider a continuous variable -order control. Again, since we are assuming that the 
number of evaluations of the derivative governs the overall efficiency, the effects of this control on 
C (in the previous section), and its possible adverse effects on A when used in conjunction with a 
variable step control, make this control ineffective from this point of view. If, however, a fixed- 
step method is used, the effect of this control on the number of corrector iterations can result in 
considerable efficiency gains. 


NUMERICAL RESULTS 

This section discusses the results obtained by applying the various local error controls dur- 
ing the integration of three selected orbit types. Appendix B formulates the actual equations of 
motion (Equation 1); a description of the computer program can be found in Reference 8. 

For the sake of simplicity, instead of an approximation for Rn (given by Equation 6) in control- 
ling the error, a bound on the local error given by 


Un = Ch 2 V k ~ 2 x 

n 


was used, where V k " 2 is the last backward difference retained in the computations (see Appendix A). 
Since p = k + 1, therefore 


I Rn | < i Un | 

Thus using Un gives the same qualitative results as using Rn, and the same quantitative results 
could be obtained by using Rn with a smaller T 1 , 

The number of derivative evaluations for each integration was obtained by multiplying the total 
number of steps taken by the average number of predictor-corrector iterations. 

The error estimates were obtained by integrating Equation 1 with P = 0, and comparing the 
results with the Kepler solution. Since the perturbation function p generally has only a small effect 
on the accumulation of error, the error estimated in this manner can be considered close to the 
actual error. 

All computations were performed on the Univac 1108 computer in double precision. The fol- 
lowing units were used: 

Unit of length = 6378.388 km, 

Unit of time = 13.447 min. 
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Table 1 gives the results obtained by integrating the orbits with formulas of orders 7, 9, 11, 
and 13. These integrations were carried out with various stepsizes, where in each case the step- 
size was increased until Un failed to satisfy the inequality 

lUnl < Tj (12) 


over the entire range of integration (chosen arbitrarily as 4000 minutes). The table gives only the 
last 2 stepsizes tested, indicating the largest stepsize that passed the criterion, and the first step- 
size failed. 


Semi- 
major 
Axis, a 


6.7 


1.15 


8.5 


Table 1 


Fixed-Order, Fixed-Step Control. 


Tj = 0.5 x 10 8 S = 0.1 x 10" 10 (predictor-corrector tolerance) 


Eccentricity 

e 


0.003 


0.075 


0.878 


Order 

Stepsize 

(min) 

No. of 
Derivative 
Evaluations 

Estimated 

Error 

7 

5.0 

798 

4 x 10' 10 


7.0* 

569 

4 x 10' 9 

9 

15.0 

262 

3 x 10 -9 


17.0* 

231 

8 x 1CT 9 

11 

24.0 

160 

<D 

X 

H-* 

O 

1 

o 


26.0* 

147 

2 x 10' 9 

13 

22.0 

173 

3 x 10' 12 


24.0* 

158 

9 x 10' 11 

7 

0.4 

9998 

8 x 10' 8 


0.5* 

7998 

1 

o 

1 — 1 

X 

9 

0.9 

4440 

2 x 10" 7 


1.0* 

3996 

5 x 1<T 7 

11 

1.2 

3327 

5 x 10' 8 


1.3* 

3070 

1 x 10' 7 

13 

1.5 

3081 

1 x 10~ 9 


1.6* 

3043 

1 x 10' 8 

7 

0.30 

13344 

1 x 10' 7 


0.35* 

11483 

X 

H 1 

o 

9 

0.40 

10005 

00 

1 

o 

T— 1 

X 


0.45* 

8902 

1 X 10' 7 

11 

0.30 

13340 

9 x 10 -11 


0.35* 

11433 

4 x 1CT 10 

13 

0.20 

19992 

4 x 1<T 12 


0.25* 

15992 

3 x 10" 12 


Integration with this stepsize failed the criterion (12) for some n. 
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The results in Table 1 suggest the following observations: 

For a given local error criterion, increasing the order does not generally imply a larger 
"allowable” stepsize. Since the propagation of error (both truncation and round-off) influences the 
"smoothness" of thd higher-order differences (and hence our local error estimate), this behavior 
is to be expected, both from the inaccuracies in the differences and from an unstable P and h com- 
bination. Also, accuracy increases with the higher orders. Reference 6 includes a table demon- 
strating this behavior (for the case e = 0). 

The stepsizes resulting from the local error criterion during the integrations correspond to 
bounding the local error at its maximum (at perigee), although the local error for most of the in- 
tegration, particularly for the case e = 0.87, was considerably smaller than the upper bound. 

Table 2 gives the results obtained by integrating the test orbits with a variable-step control. 

In particular, the stepsize was varied during each integration, forcing Un to satisfy 

Tl < |Un| < T 2 (13) 

for all n. The stepsize- modification techniques used were halving-doubling and the optimum-step 
computation where the "allowable" local error is designated by a-; i.e., the step computation is 
given by 


>V t = [orh p+ 2 /Un ] 1/p+2 . 

Except for those cases indicated by an asterisk, all integrations were performed with an initial 
stepsize (chosen arbitrarily) of 0.42 min = 1/2 5 internal units. The sixth column indicates the 
stepsize or stepsize range that occurred during the multistep integration as a result of the local 
error control. 

Comparing the corresponding results in Tables 1 and 2 leads to the following observations: 

For orbits of low eccentricity, (0.003, 0.075), the stepsizes selected in Table 2 as a result of 
the variable-step control are comparable with the "optimum" stepsizes found by trial and error in 
Table 1. 

Table 2 indicates by asterisks those cases in which the initial step was selected as a result of 
foreknowledge of the "optimum" step (the initial step used was exactly half the stepsize indicated). 
These cases demonstrate the dependence of the efficiency on the initial choice of stepsize in the 
halving-doubling type of step modification. 

For orbits of high eccentricity (0.87), the larger stepsizes effected significant gains in efficiency 
Stepsize modification causes no loss in efficiency, as shown by the results concerning the degree of 
control in Table 2. The apparent loss in accuracy, (especially for the higher orders) was expected 

10 
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Table 2 


Fixed-Order, Variable-Step Control. 


Semi- 
major 
Axis, a 


6.7 


1.15 


8.5 


Eccentricity 

e 

0.003 


0.075 


0.87 


T x = 0.5 x 10“ 8 T 2 = 0.5 x HT 13 8 = 0.1 x 10~ 10 


Order 

P 

Type of 
Step 

Modification 

O' 

Step size 
(Range) 
(min) 

No. of 
Derivative 
Evaluations 

Estimated 

Error 

7 

H/D 


1.7 

2377 

2 x 10 “ 13 


OPT 

H- 1 

X 

h- 2 

O 

1 

O 

5.4 

744 

6 x 10“ 10 


OPT 

1 x 10" 11 

3.0 

1320 

1 x 10“ 11 

9 

H/D 


6.7 

590 

2 x 10 “ 12 


H/D* 


14.0 

281 

1 xlO" 9 


OPT 

1 x 10“ 10 

12.4 

322 

5 x 10 “ 10 


OPT 

1 x 10' 11 

8.4 

471 

1 x 10“ n 

11 

H/D 


13.4 

291 

2 x 10“ 12 


H/D* 


23.0 

167 

6 x 10 “ 10 


OPT 

1 x 10“ 10 

20.1 

196 

1 x 10 ” 10 


OPT 

1 x 10“ n 

15.1 

260 

7 x 10“ 12 

13 

H/D 


13.4 - 26.9 

174 

6 xl0“ 9 


H/D* 


23.0 

165 

5 x 10“ 12 


OPT 

1 x 10' 10 

19.8 - 30.0 

155 

8 x 10 “ 9 


OPT 

1 X 10' 11 

13.9 - 26.8 

159 

7 x 10“ 9 

7 

H/D 


0.42 

9516 

1 x 10 “ 7 


OPT 

1 x 10" 10 

0.42 

9516 

1 x 10 7 


OPT 

1 x NT 11 

0.42 

9516 

1 x 10 ” 7 

9 

H/D 


0.42 

9514 

2 x 10 “ 10 


OPT 

i x io -10 

0.61 

6514 

6 x 10 “ 9 


OPT 

1 x 10' 11 

0.42 

9490 

2 x 10 “ 10 

11 

H/D 


0.84 

4781 

1 x 10 “ 9 


OPT 

1 X 10“ 10 

1.04 

3849 

8 xl0“ 9 


OPT 

1 x 10” 1 1 

0.78 

5131 

4 xlO“ 10 

13 

H/D 


0.84-1.68 

4257 

4 xl0“ 7 


OPT 

1 xlO' 10 

1.20 

3314 

8 x 10 “ 10 


OPT 

i x nr 11 

0.96 

4171 

5 x 10 “ 10 

7 

H/D 


0.2 -3.4 

3180 

7 xl0“ 9 


OPT 

1 x 10* 10 

0.2 -9.5 

2415 

2 x 10“ 8 


OPT 

1 X IO' 11 

0.1 - 6.9 

3131 

3 x 10“ 9 

9 

H/D 


0.2-6. 7 

1374 

7 x 10“ 8 


OPT 

X X 10“ 10 

0.2-19.0 

1137 

7 x 10 “ 8 


OPT 

1 x 10 “ 1 1 

0.2-14.6 

1331 

1 x 10“ 9 

11 

H/D 


0.2-13.4 

875 

6 x 10 “ s 


OPT 

1 x 10“ 10 

0.2-29.1 

788 

2 x 10 “ 8 


OPT 

1 x 10“ 1 1 

0.2-14.0 

907 

5 x 10 “ 9 

13 

H/D 


0.2-26.9 

710 

1 x 10 “ 7 


OPT 

1 x 10 “ 10 

0.2-22.4 

661 

3 x 10 “ 8 


OPT 

1 x 10“ 1 1 

0.1-25.6 

775 

1 x 10 “ 8 


H/D - Halving-doubling 

OPT “ Optimum-Step Computation 
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since the error control bounded the error from below as well as from above, whereas for the cor- 
responding results in Table 1 the local error was insignificant for most of the integration. 

In all cases, the optimum-step computation-type interval modification yielded fewer derivative 
evaluations than halving-doubling, with little or no loss of accuracy. The parameter magnitude af- 
fected the integration as expected: the small value increased the accuracy at the cost of smaller 
stepsizes. 

In all cases, automatic selection of a stepsize to satisfy the local error criterion for the given 
order is clearly an advantage. Also, the "best" order (i.e., the one yielding the fewest derivative 
evaluations) to use with the variable-step control may not be the same as the best order to use for 
a fixed-order, fixed- step integration. 

Table 3 gives the results obtained by integrating the test orbits with a variable-order, variable- 
step control. The order was allowed to vary between the limits L x and l 2 before any step modifica- 
tion was effected. In all cases, the "smallest" order satisfying Equation 12 in any given stepsize 
was used. All integrations were performed with an initial stepsize of 0.42 mins, and an initial order 
of Lj + (l 2 -L^/2, forcing the initial order p to be in the interval (l 15 l 2 ). Columns 5 and 6 indicate 
the stepsize and order ranges that occurred during integration as a result of the error control. 

Table 3 

Variable -Order, Variable-Step Control. 



Tj = 0.5 

o 

i — i 

X 

-8 

o 

i — i 

X 

LO 

o 

tl 

CN 

H 

8 = 0.1 x 

0 

1 

O 

7 = 0.-1 X 10 " 9 


Semi- 
major 
Axis, a 

Eccentricity 

e 


" L 2 

Type of 
Step 

Modification 

Stepsize 

(Range) 

(min) 

Order 

(Range) 

No. of 
Derivative 
Evaluations 

E stimated 
Error 

6.7 

0.003 

7 

- 13 

H/D 

3.3 

8 

1184 

9 x 10" 14 





OPT 

8.8 

10 

454 

5 x 10" 13 



9 

- 13 

H/D 

6.7 

9 

588 

2 x 10" 12 





OPT 

18.4 

11 

217 

5 x 10" 11 



11 

- 15 

H/D 

13.4 

11 

289 

2 x 10" 12 





OPT 

24.6 

11 - 15 

175 

2 x 10" 9 

1.15 

0.075 

7 

- 13 | 

H/D 

0.42 

8-10 

9513 

9 x 10“ 10 





OPT 

0.84 

10 

4767 

5 x 10 -9 



9 

- 13 

H/D 

0.42 - 0.84 

9-10 

4818 

6 x 10" 9 


i 

i 



OPT 

1.04 

11 

3849 

8 x 10" 9 



11 

- 15 

H/D 

0.84 

12 - 13 

4751 

1 x 10“ 10 





OPT 

1.2 

13 

3314 

8 x 10' 10 

8.5 

0.87 

7 

- 13 

H/D 

0.2 - 3.3 

7-13 

2744 

4 x 10" 9 





OPT 

0.2 - 10.8 

7-13 

2269 

1 x 10" 8 



9 • 

- 13 

H/D 

0.2 - 13.4 

9-13 

1002 

6 x 10" 8 





OPT 

0.2 - 13.2 

9-13 

1067 

3 x 10" 8 



11 

- 15 

H/D 

0.1 - 26.8 

11 - 15 

665 

3 x 10" 7 





OPT 

0.1 - 20.0 

11 - 15 

742 

2 x 10" 7 
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Comparing Table 3 with Tables 1 and 2 suggests the following observations: 

For orbits of low eccentricity, in Table 3, the stepsizes attained for the various orders selected 
were comparable to the stepsizes obtained for the corresponding orders in Table 2. 

In all cases, (especially e = 0.87), allowing the order to vary before changing the stepsize re- 
sulted in a smaller "mean” stepsize and thus a larger number of derivative evaluations; so that if 
the "best" order to use with a variable step control is known, this control gives a more efficient 
integration than the variable-order, variable-step control. . 

In all cases, automatic selection of a stepsize and order to satisfy the given local error crite- 
rion is clearly an advantage. 

Table 4 shows the effects of a variable-order control alone on the average number of predictor- 
corrector iterations (and hence the total number of derivative evaluations) and on the total error. 

The order or order range obtained as a result of the control is indicated. 

Table 4 

Variable -Order, Fixed-Step and Fixed-Order, Fixed-Step Control. 


T t = 0.5 x 10 -8 T 2 - 0.5 x 10~ 13 8 - 0.1 x 10 -10 


Semi- 
major 
Axis, a 

Eccentricity 

e 

Mode 

Stepsize 

Order 

No. of 
Steps 

Average Number 
of p-c 
Iterations 

Estimated 

Error 

6.7 

0.003 

Vary order 

25.0 

11 

154 

1.00 

1 X 10* 9 



Fixed order 

25.0 

7 

158 

1.80 

4 x 10* 6 



Fixed order 

25.0 

9 

156 

1.00 

2 x 10* 7 

1.15 

0.075 

j Vary order 

! 1.2 

11 

3327 

1.00 

5 x 10" 8 



Fixed order 

1.2 

7 

3331 

1.94 

1 x 10* 5 



Fixed order 

1.2 

9 

3329 

1.28 

2 x 10" 6 

8.5 

0.87 

Vary order 

0.5 

5-10 

8000 

1.00 

2 x 10* 8 



Fixed order 

0.5 

7 

7998 

1.01 

\D 

1 

o 

I— I 

X 

LO 


These results show the advantages of automatic selection of order, from the predictor-corrector 
point of view and from the accumulated error that results from local error control. 

Finally, consider the effects of varying the "allowable range" (T lt t 2 ) in the local error con- 
trol. In general, the smaller the interval (t 1? T 2 ) is made, the greater is the frequency of interval 
and/or order modification. It would seem (particularly in the case of a variable-step control), that 
efficiency gains due to the larger stepsizes due to the control could be offset by the cost in chang- 
ing the stepsize— should such changes occur too frequently during integration. 

Table 5 gives the results of integrating one of the test orbits with a fixed T x and various values 
for T 2 (approaching T 1 as a limiting value). A variable-step control (halving-doubling) was used, 
and the variations in the stepsize, along with the number of step changes and computation time (in 
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Table 5 


Variable-Step Control — Range Modification 


Tj = 0.5 x 10~ 8 S = 0.1 x 10~ 10 


Semi- 
major 
Axis, a 

Eccentricity 

e 

T 2 

Stepsize 

(Range) 

(min) 

No. of 
Derivative 
Evaluations 

8.5 

0.87 

5 

x 10' 16 • 

0.42 - 6.7 

1810 



5 

x 10‘ 14 

0.42 - 6.7 

1309 



5 

x 1CT 12 

0.42 - 13.4 

886 



1 

x 10' 11 

0.42 - 13.4 

849 



2 

x 1CT 11 

0.42 - 13.4 

832 



2.5 

X 

O 

t 

0.42 - 26.9 

856 



3 

* 

1 

o 

I— 1 

X 

0.42 - 26.9 

892 


Estimated 

Error 


7 x 10" 8 
7 x 10~ 8 
6 x 10' 8 
9 x 10" 8 
1 x 10" 7 
lxlO' 7 
1 x 10" 7 


Computation 

Time 

(min) 

0.094 

0.073 

0.059 

0.058 

0.260 

0.810 

2 . 10 * 


minutes) are presented for each integration. An asterisk indicates that integration in which a step- 
size modification occurred approximately at each step, rendering the error control completely in- 
effective, since the entire computation time was governed by the backpoint computation. 


These results show that as T 2 approaches there are efficiency gains at first, because of the 
larger stepsizes, but, as the range ( T 1 , T 2 ) becomes smaller, these gains are offset by the cost in 
changing the step. Finally, as expected, the diminishing range causes the control to be completely 
ineffective. Note, however, that the equations of motion (see Appendix B) are simpler than what 
may be used in actual practice, and in a more realistic situation, a smaller interval (t 1? t 2 ) may 
be possible before the computing time is completely governed by the frequency of step modification. 


It may be admitted that only a particular model of the perturbation function P has been used. 
However, variations in this function (such as the inclusion of higher-order gravitational effects) 
should not affect the general qualitative behavior of the local error. 


CONCLUSIONS 

The problem of achieving an efficient and accurate orbital integration by a multistep process, 
using the concept of controlling the local error during integration, has been studied. It has been 
shown that during the integration, the parameters p and h can be used to control the local error in 
such a way that the efficiency of the process is improved with no significant loss in accuracy. 

In particular, if a sufficiently large range (t 1? T 2 ) in the local error is allowed, and an order 
p is given, a variable- step control can automatically yield a good approximation of the optimal 
initial stepsize (with respect to the given order) and can significantly improve the efficiency of the 
process by varying this step during integration. Moreover, if a good approximation of the ’’best” 
order to use with a variable -step, fixed-order control is known, this control effects the most 
efficient integration, as compared with the other controls considered. On the other hand, a 
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variable-order, variable-step control automatically determines both the order and step required 
to satisfy the given local error criterion, and also gives a reasonably efficient integration. Finally, 
even a variable-order, fixed-step control, although it does not minimize the number of integration 
steps, can give a more efficient integration than no control at all, by minimizing the number of 
predictor-corrector iterations. 

We have considered only a ’local" optimization problem, in the sense that optimal stepsize is 
defined on a step-by-step basis as being the largest stepsize satisfying a given local error criterion 
at any given point. A more significant consideration would be optimal stepsize (and/or order) dis- 
tribution over the entire range of integration, defined on the basis of a criterion imposed on the ac- 
cumulated truncation error. D. Morrison, under some restrictive conditions (among others, limi- 
tation to a single differential equation), considers the problem of optimizing the mesh distribution 
(Reference 9). A basic difficulty in applying such techniques is the requirement that a "memory 11 
of equally spaced points be available at each step during the integration. One solution worth con- 
sidering is the integration of divided-difference interpolation polynomials (which do not require 
equally spaced points), for use in numerical integration. 


Goddard Space Flight Center 
National Aeronautics and Space Administration 
Greenbelt, Maryland, January 22, 1968 
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Appendix A 

Derivation of the "Summed” Form of Integration Formulas 


By integrating Newton's backward difference interpolation polynomial, the following multistep 
formulas for the numerical integration of Equation 1 can be derived,* (including formulas for the 
velocity in the event that p contains x): 

Equations for x x 


Predictor: x - 2x m + x . 

1 m mi 

Corrector: x - 2x + x m . 

m+ 1 rn m— 1 

Equations for x - x 

Predictor: x - x 

m'r 1 m 

Corrector: x - x m 

m+ 1 m 


+ CT- V 3 X + * ' * + Cr V k X 1 , 
n 3 m k m J 

(Al) 

v2 ^i + --- + < vk ^ m+ i] • 

(A2) 

' k vk *J ■ 

(A3) 

■ + ?k VkW m+l] • 

(A4) 


where the coefficients cr, o% y, y* are given by the following recurrence relationships: letting 


a o = 1 


♦Reference X, pp. 192-195 and 291-293. 


m 

U 

j=l 


y* = 1 


m oh 

~ ’ zn i+l 


-L t 

j=l 


j + 2 


2h j+i , 

‘ L i + 2 ^ 

j=i 
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^m 



j=l 



j=l 


m = 1 , 2 , 3 , * * • , k , 


and 

V *m = *m ' *.-l ■ 

V 2 X = V (Vx ) = v(x "X ) = X - 2x + X , 

V k X = vfv * 1 - 1 X ) . 

It has been established (References 1 and 10) that an algebraic equivalent of Equations A1 
through A4, known as the ’’summed” form of the integration formulas, considerably reduces the 
propagation of round-off error. Formally, one can obtain this equation modification by applying 
the inverse difference operators V’ 1 , V" 2 defined by V " 1 V = I , V " 2 V 2 = i , (i the identity), to both 
sides of these equations. In particular, applying V' 2 to both sides of Equations A1 and A2 gives 


Predictor : x +1 h 2 V ~ 2 x + a V“ 1 x + cr x + a Vx + * * • + cr V k 2 x 1 , 

m+l L 0 ml m 2 m 3 m k mJ 


Corrector: x m+1 - h 2 [cr Q * V " 2 x m+1 +cr* V ’ 1 x m+1 +a* x^ + a* Vx m+l + • • • +a k * V k " 2 x m+ J 


Defining ! s , n S by 


J S = V _I ( : S ) = V' 2 i 

m \ m / ir 


gives 


Predictor: x ,, 

m+l 


\a n ng 

10 m 


+ a. 


+ cr x 

i 2 m 


hcr t V k_2 x 1 , 

k m J 


(A5) 
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and since 


*s 


*m+l 


+ *s 


and 


n s 


n S m 


: S + x 


therefore 


Corrector: x m+ , - h 2 [o- 0 * n S m + (a* + cr* ) I S m + (c 


+ a. 


+ Cr 2*) *«+l + cr * 


Vx 


+ cr,*V k * 2 : 


Ui](A6) 


Likewise for Equations A3 and A4, applying the operator V 1 to both sides gives 


Predictor: i m+1 - h [y Q V 1 x m + y 1 x m + y 2 Vx m + • • • + y k V k ' 1 x m ] , 


Corrector: i m+1 = h [y 0 * V ' 1 x m+1 + y* x m+1 + y 2 * Vx m+1 + • • • + y* V^i x m+1 ] 


Using the definition of T s as above gives 


Predictor: x m+1 = h [y 0 ^ + y, x m + y 2 Vx m + • • • + y k V k "‘ x m ] , 


(A7) 


Corrector: i m+1 = h [y 0 * I S_ n + (y 0 * +y*) x m+1 + y 2 * ^ m+1 + • • • 


1 X m + l] • (A8) 


These are the ” summed” forms of the integration formulas. (For details concerning the computa- 
tional usage of these formulas, see Reference 8.) 

Note, again, that Equations Al, A2 and A5, A6 are algebraically equivalent, so that for any 
fixed k, any solution of Al, A2 is a solution of A5, A6, and conversely. In particular, the local 
truncation errors associated with these formulas are the same. A similar equivalence exists for 
Equations A3, A4, and A7, A8. 
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Appendix B 

Equations of Motion 


The equations of motion used to obtain the numerical results are given by 


if* + f 

t?3 M 


+ F 0 


where F 1 + F 2 = P (see Equation 1), and x = (x, y, z) , R = I x| = (x 2 + y 2 + z 2 ) 1/2 , ^ is a constant, F 1 
is the perturbation due to the non- sphericity of the earth, and F 2 is the pertu^ 1 ation due to drag. 

F i = ( F ix- F iy- F iz) is given by 

F lx = g [jR 2 (5s- 1) + Hz(7s- 3)+ § (-63S 2 +42s -3)] , 

F i y = ^[jR 2 (5s-l) + Hz(7s-3)+£(-63s 2 +42s-3)] , 

F t, " g[jR 2 (5s-l) + |-(-63s 2 +70s-15)] + ^ (35s 2 -30s+3) , 

where s = (z/R) 2 ; and j, H, and K are the second, third, and fourth harmonics of the earth T s po- 
tential field, respectively. f 2 is of the form 

C D A 

arplvj v r , 

where V r is the relative velocity of the satellite; p is the atmospheric density at the satellite po- 
sition; and A, M and C D are the cross-sectional area, mass, and drag coefficient of the satellite re- 
spectively. The actual numerical values of these constants that were used to obtain the results can 
be found in Reference 8. 
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